arXiv:astro-ph/9701140vl 20 Jan 1997 


Highly Viscous Accretion Disks with Advection 


I.V.Artemova^, 

Theoretical Astrophysics Center, Juliane Maries Vej 30, DK-2100, Copenhagen 0, Denmark 

G.S. Bisnovatyi-Kogan^, 

Space Research Institnte, Profsoyuznaya 84/32, 117810 Moscow, Rnssia 

G. Bjornsson^, 

Science Institnte, Dnnhagi 3, University of Iceland, IS-107 Reykjavik, Iceland 

I.D. Novikov'^ 

Theoretical Astrophysics Genter, Juliane Maries Vej 30, DK-2100, Gopenhagen 0, Denmark 
University Observatory, Juliane Maries Vej 30, DK-2100, Gopenhagen 0, Denmark 
NORDITA, Blegdamsvej 17, DK-2100 Gopenhagen 0, Denmark 
Astro Space Genter of P.N. Lebedev Physical Institute, Profsoyuznaya 84/32, 117810 

Moscow, Russia 

Received_; accepted_ 


^e-mail:julia@nordita.dk 
^e-mail; Gkogan@mx.iki.rssi.ru 
^e-mail;gulli@raunvis.hi.is 
“^e-mail; novikov@nordita.dk 





2 


ABSTRACT 

We consider the effects of advection and radial gradients of pressure and 
radial drift velocity on the structure of optically thick accretion disks. We 
concentrate our efforts on highly viscous disks, a = 1.0, with large accretion 
rates. Contrary to disk models neglecting advection, we hnd that continuous 
solutions extending from the outer disk regions to the inner edge exist for 
all accretion rates we have considered. We show that the sonic point moves 
outward with increasing accretion rate, and that in the innermost disk region 
advection acts as a heating process that may even dominate over dissipative 
heating. Despite the importance of advection on it’s structure, the disk remains 
geometrically thin. 


Subject headings: accretion, accretion disks - black hole physics 
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1. Introduction 

The “standard accretion disk model” of Shakura (1972) and Shakura & Sunyaev (1973), 
that has been widely used to model accretion flows around black holes, is based on a 
number of simplifying assumptions. In particular, the flow is assumed to be geometrically 
thin and with a Keplerian angular velocity distribution. This assumption allows gradient 
terms in the differential equations describing the flow to be neglected, reducing them to 
a set of algebraic equations, and thereby hxes the angular momentum distribution of the 
flow. For low accretion rates, M, this assumption is generally considered to be reasonable. 

Since the end of the seventies, however, it has been realized that for high accretion 
rates, advection of energy with the flow can crucially modify the properties of the innermost 
parts of accretion disks around black holes. A deviation from a Keplerian rotation may 
result. 

Initial attempts to solve the more general disk problem only included advection of 
energy and the radial gradient of pressure in models with small values of the viscosity 
parameter, a = 10“^ (Paczyhski & Bisnovatyi-Kogan 1981), and it was shown that 
including radial velocity in the radial momentum equation would not change principally 
the results for such a small a (Muchotrzeb & Paczyhski 1982). Liang & Thomson (1980) 
emphasised the importance of the transonic nature of the radial drift velocity, and the 
influence of viscosity on the transonic accretion disk solutions was noted by Muchotrzeb 
(1983), who claimed that such solutions only existed for viscosity parameters smaller than 
a* ~ 0.02 — 0.05. Matsumoto et ah (1984), then showed that solutions with a > a* do 
in fact exist, but the nature of the singular point, where the radial velocity equals the 
sound velocity, is changing from a saddle to a nodal type and the position of this point 
is shifted substantially outwards in the disk. Matsumoto et al. (1984) also demonstrated 
the non-uniqueness of the solutions with a nodal type critical point for given Keplerian 



4 


boundary conditions at the outer boundary of the disk (see also Muchotrzeb-Czerny 1986). 
Extensive investigation of accretion disk models with advection for a wide range of the disk 
parameters, M and a, was conducted by Abramowicz et ah (1988), with special emphasis 
on low a. Misra & Melia (1996) considered optically thin two-temperature disk models 
and treated advection in the framework of the Keplerian disk model, but hxed the proton 
temperature somewhat arbitrarily at the outer boundary. Chakrabarti (1996) solved the 
advection problem containing shock waves near the innermost disk region, considering 
accretion through saddle points. Numerical solutions of accretion disks with advection have 
been obtained by Chen and Taam (1993) for the optically thick case with a = 0.1, and by 
Chen et ah (1996), for the optically thin case (see also Narayan 1996). A simplihed account 
of advection has recently been attempted, either treating it like an additional algebraic 
term assuming a constant radial gradient of entropy (Abramowicz et ah 1995; Chen et ah 
1995; Chen 1995), or using the condition of self-similarity (Narayan & Yi 1994). 

Over the last few years it has become clear, that neglecting the advective heat transport 
at high M leads to qualitatively wrong conclusions abont the topology of the family of 
solutions of the disk structnre eqnations (see for example Abramowicz et ah 1995; Chen 
et ah 1995; Artemova et ah 1996). The disk structure equations withont advection give 
rise two branches of solutions: optically thick and optically thin, which do not intersect if 
M < Mfe Ri (0.6 — 0.9)MEdd for a = 1 and Mbh = IO^Mq, where iifEdd is the Eddington 
accretion rate (Artemova et al. 1996). For larger accretion rates there are no solntions 
of these equations extending continuously from large to small radii, and with Keplerian 
bonndary conditions at the onter bonndary of the disk (see also Liang & Wandel 1991; 
Wandel & Liang 1991; Luo & Liang 1994). It was argued by Artemova et al. (1996), that 
for accretion rates larger than Mf,, advection becomes critically important and would allow 
solutions extending all the way to the inner disk edge also to exist for M > M^. 
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The goal of the present paper is to construct explicitly accretion disk models for 
high M and large a taking advective heat transport self-consistently into account. We 
also include radial gradients of pressure and radial drift velocity and we allow for the 
non-Keplerian character of the circular velocity. Furthermore, we use the geometrically thin 
disk approximation because, as will be seen in our solutions, the relative thickness of the 
disk is substantially less than unity. We show that solutions extending from large radii to 
the inner edge of the disk can be constructed even for accretion rates considerably larger 
than Mfe. We hnd that advection is very important in the innermost disk region, although 
the flow does not deviate strongly from Keplerian down to the region where the radial 
inflow velocity approaches the local sound speed. 

In §2 we introduce our model and describe our solution methods, while in §3 we discuss 
our results. 


2. The Model and the Method of Solution 

In this paper we will only consider optically thick solutions to the disk equations. When 
advective cooling is important we assume that it can be sufficiently well modelled by adding 
it self-consistently to other cooling mechanisms in a geometrically thin disk. 

We use from now on geometric units with G = 1, c = 1, use r as the radial coordinate 
scaled to Vg = M, and scale all velocities to c. We work with the pseudo-Newtonian 
potential proposed by Paczyhski and Wiita (1980), <h = —M/(r — 2), that provides an 
accurate, yet simple approximation to the Schwarzschild geometry. We normalise the 
accretion rate as m = M/M-Edd, where M^dd = L^dd = 4:nMmp/aT, in our units. 

We use the same equations and ingredients in our models as in Artemova et al. (1996), 
except for changes required by the Paczyhski-Wiita potential and the differential terms 
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in the energy equation and the radial momentum equation. The following equations are 
therefore modihed: 


1. Conservation of angular momentum for a steady-state accretion in the a-disk model, 
is written as 


rh (ugfi) 


dlnfl 


dlnr 


-1 


/ = 


(Tt 


nir 


haP, 


( 1 ) 


where, fl, is the angular velocity, the factor / = 1 — /in//, where, / = is the specihc 
angular momentum and /in is the value of / lost from the disk at the innermost edge and 
swallowed by the black hole. The half thickness of the disk is denoted by h, and P is the 
total pressure in the equatorial plane of the disk. 


2. The energy equation has the form 


Q+ Qa,dv T Ql oc: 


( 2 ) 


where Qioc is in general the total rate of all local cooling processes (see Artemova et ah 
1996), and the viscous heating rate per unit area, Q+, is given by the formula (see e.g. 
Bisnovatyi-Kogan 1989; Frank, King & Raine, 1992) 
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The advective cooling rate can be written in the form (see e.g. Chen & Taam 1993): 
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( 4 ) 


27rr dr 27rr 

where T is the temperature and S is the specihc entropy. Here, E is the energy per unit 
mass of the gas and v = 1/p, where p is the matter density. With Qadv of the form given in 
equation (^), the energy balance becomes a differential equation. 

3. The momentum equation in the radial direction takes into account pressure and 
radial velocity gradients and is written in the form: 


\ d dll 

-^ = {n^-nl)r-vr^, 

p dr dr 


( 5 ) 
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where, f^K = \J{d^/dr)/r, is the Keplerian angular velocity in the Paczyhski-Wiita 
potential and Vr is the radial drift velocity. Neglecting the gradient terms in equation (^, 
as is done in the standard model, hxes 17 = 17 k- 


From equation (|^) and mass conservation one gets: 

dlnvr 05(1 + {dlnh)/{dlnr)) + r^(17^ — 17^) 
d In r — al 


where is dehned as (Muchotrzeb & Paczyhski 1981): 


a 


2 

S 


f(W\ ^ 

\dr ) \dr j 


( 6 ) 

(7) 


Note, that a* is not a physical sound velocity but rather a formal quantity. The vanishing 
of both the numerator and the denominator in equation (|^) at the same value of r, the 
’’sonic point”, provides the regularity condition required for a ’’transonic solution” of the 
flow structure. 


To solve the differential equations (|^) and (||) we adopt the following boundary 
conditions: At large radii, r S> 100, the solution must coincide with the standard Keplerian 
disk solution obtained neglecting advection. In addition, the parameter Zjn in equations (|^ 
and (]^), is an eigenvalue of the problem which is adjusted in such a way that the solution 
satishes the regularity condition at the “sonic point”. 

We solved this system of equations numerically by the method of subsequent iterations 
with hxed rh and a. Starting from the ’’standard disk” solution as the initial trial solution 
with a specihc value of in the function, / = 1 — kn/l, we then varied to obtain a 
self-consistent solution. Typically the method converges to a solution after three to four 
iterations. 


In practice, we varied /in in some interval and determined the positions of the points 
r^v and r£, where the numerator and the denominator of equation vanish, respectively. 





We then considered the dependence of the difference — r^r) on /in and determined /in for 
which the difference (r^ — r^r) is equal to zero. The corresponding /in is an eigenvalue of the 
problem. Examples of the dependences {ro — r^) on /in are given in Figure 1 for m = 10.0 
(accretion rate less than rhb), and rh = 28.0 (accretion rate substantially greater than rhb). 

For m > TTib, our method is very sensitive to the choice of the initial trial solution in 
the innermost disk region. Using a “standard disk” solution down to r = 6 is not possible, 
as for these large values of m there is no solution around r ^ 13, and the method cannot 
bridge that gap to hnd a “transonic” solution. For our initial trial solution, we therefore 
chose a value of /in that allowed us to generate the trial solution down to small radii, and 
then iterated as described above. 

As is seen in Fig. 1, there are three values of /in for each m, where (r^) — r^v) = 0, and 
some range of /in where {td — tn) is very close to zero. Most likely that range corresponds 
to the nodal type of a sonic point at large a, as obtained by Matsumoto et al. (1984) and 
others under some simplihcations. In this case the condition of regularity at the sonic point 
does not specify /^ uniquely. From our numerical method we are unable to determine if 
any /^ in the range of /in’s, where — vn is close to zero, provides an acceptable solution. 
Complete analysis of the character of the critical points needs a different approach and will 
be performed elsewhere (but see next section). 

Our method allows us to construct a self-consistent solution to the system of equations 
from very large radii, r > 100, and down to the innermost regions of the disk. The solution 
passes through a ’’sonic point” and continues closer towards the black hole. But, we cannot 
construct the parts of the solutions in the very vicinity of the black hole where the angular 
velocity is very far from Keplerian and / is almost constant. However, in all cases do we 
extend the solutions down to the value of r at which Vr becomes equal to the local adiabatic 
sound velocity. These radii are in general closer to the black hole than the location of the 



sonic point”. 


3. Results and Discussion 

In Table 1 we summarise the parameters of the models for which — tat = 0, according 
to our computations. For each fixed m, the properties of the two (or three) self-consistent 
solutions are similar and differ only quantitatively. In all cases discussed below do we take 
Mbh = 10 ®Mq and a = 1. 

We will now compare the solutions with and without advection. In the “standard 
model”, for accretion rates m < mh = 14.315, there always exist solutions that extend 
continuously from large to small radii. When m > rhb = 14.315 there are no solutions in a 
range of radii around r ^ 13, and therefore no continuous solutions extending from large 
radii to the innermost disk edge (see detailed discussion by Artemova et ah 1996, where 
however, the Newtonian potential was used, resulting in rhb = 9.4). 

In Figure 2a we plot the disk surface density S = 2ph, in an optically thick disk as 
a function of radius, r, in a model with m = 10. The lowest curve is the solution of the 
standard model, the upper ones are solutions number 2 and 3 in Table 1. Note that the 
solutions including advection all terminate at radii considerably greater than r = 6 (inner 
edge of the disk in the standard model). 

In Fig. 2b we plot similarly the solutions for m = 15. In the standard model, no 
solution exists in the region around r si 13, but when advection is included, the structure 
of the solutions is completely different. Models 6 and 7 in Table 1 are shown. 

For m < 13, including the gradient terms gives rather small corrections to the standard 
disk model, see Figure 2a. When rh > 13 advection becomes essential and for m > hib it 
changes the picture qualitatively. When rh > rhb solutions do exist extending continuously 
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from large radii to the innermost disk region where the solution passes through a ’’sonic 
point” (compare Figs 2a and 2b, see also Fig 4b below). As mentioned above (see the 
end of Section 2), we can extend our models only down to the region where the radial 
velocity becomes equal to the local adiabatic sound velocity. We are unable to calculate the 
properties of the flow for smaller radii. Only more detailed analysis of this region (using 
other methods) allows one to determine the smoothness of the flow down to the event 
horizon of a black hole or verify the presence or absence of shocks in the region. 

In Figure 3 we plot a family of optically thick solutions for different m, clearly 
demonstrating that the solutions to the complete system of disk structure equations 
including advection and radial gradients have quite different properties at high m compared 
to the solutions of the standard disk model. 

In Figure 4a we plot Qadv/Q+ as a function of radius for rh = 10 and m = 28, that 
bracket the cases we have studied. Outside the radius where the entropy gradient is zero 
(and therefore Qadv = 0, recall eq. [^), advection provides an additional cooling, that is 
however, never substantial in our models. On the other hand, inside that radius, advection 
acts as a heating process that easily dominates over the dissipation rate that decreases 
rapidly near the inner edge of the disk (as / —>• 0, see eq. 0). Panels 4b and 4c show 
the corresponding Mach numbers and h/r-ratio, respectively. Note that although the flow 
becomes transonic in the inner region, the disk can still be considered geometrically thin. 

In our calculations, the non-uniqueness of solutions at large a > a*, passing through 
the critical point (Matsumoto et ah 1984; Muchotrzeb-Czerny 1986), is preserved. It is sill 
not clear, if this non-uniqueness is a realistic physical fact which explanation may be highly 
problematic (see for example Kato et ah 1988, where the authors argue that the fact that 
the transonic point is a nodal type critical point is equivalent to an instability condition), 
or is a result of restrictive precision of our numerical solutions. Two possible approaches 
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to clarify the situation can be suggested. In the first one, we could obtain an asymptotic 
solution of the disk equations near the gravitational radius and try to match it with the 
numerical solution going from the nodal point towards the inside. The second approach 
could be hnding stationary solution by solving equations of non-stationary accretion with 
the appropriate boundary conditions. Both approaches need substantial numerical work, 
that we plan to undertake in the future. 

We would like to thank M. Abramowicz for very constructive and friendly criticism 
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FIGURE CAPTIONS 

Fig. 1. Difference (r^ — tn) as a function of Zjn, the angular momentum swallowed by the 
black hole. The 3 zero-points correspond to the solutions satisfying the regularity condition 
(eq. P), and the corresponding /jn are the eigenvalues of the problem, (a) Accretion rate of 
rh = 10 and (b) m = 28. 


Fig. 2. Disk surface density, S = 2ph, as a function of radius, comparing solutions 
with and without advection. (a) Here, rh = 10. The solid curve is the standard solution 
without advection. The dotted curve has = 3.782 and the dashed curve has kn = 4.025 
(models 2 and 3 in Table 1). (b) The case rh = 15. Again the solid curve is the standard 
model without advection. Notice the ’no solution’ region around r 13. The dotted curve 
has /in = 4.125 and the dashed curve has /in = 4.513 (models 6 and 7 in Table 1). 


Fig. 3. Surface density, S, as a function of radius for a = 1.0 and different accretion 
rates. The solid curve is for rh = 1.0, the long dashed, short dashed, dotted and dash-dotted 
curves have rh = 10,15,19 and 28, respectively (models 2, 6, 10 and 13, respectively). 

Fig. 4. Comparing the solutions for rh = 10 (dashed curve) and rh = 28 (solid curve) 
with a = 1.0, that bracket most of the cases we have studied, (a) Ratio of advective rate 
(eq. [0) to viscous heating rate (eq. 0). The advective rate equals zero when the entropy 
gradient is zero (at r 18 and r 40, for rh = 10 and rh = 28, respectively). Outside 
those radii, advection provides rather small additional cooling in both cases. Inside these 
radii advection acts as a strong source of heating, (b) Mach number and (c) ratio /r/r for 
the same cases as in panel (a). The styles of the curves are the same as in panel (a). 


Table 1. Models 


N 

m 

^in 

Ts 

1 

10 

3.570 

9.90 

2 

10 

3.782 

12.02 

3 

10 

4.025 

13.95 

4 

13 

4.010 

15.0 

5 

13 

4.170 

16.4 

6 

15 

4.125 

16.6 

7 

15 

4.513 

19.9 

8 

17 

4.251 

18.8 

9 

17 

4.850 

23.9 

10 

19 

4.416 

18.6 

11 

19 

5.088 

25.8 

12 

28 

4.305 

22.1 

13 

28 

5.169 

30.6 

14 

28 

5.409 

33.0 



























